load('storeA.mat','AA');

%Count for end lockwise(1) or anti-lockwise(2)
Count=zeros(2,1);

%Number of experiement
ExpNum=10000;

%set the start position
%The starting time is 26hr
StartPosTemp=(find(AA(:,2)>26));
StartPos=StartPosTemp(1);

TimeDisP=zeros(ExpNum,1);
TimeDisN=zeros(ExpNum,1);

for Single=1:ExpNum
    if mod(Single,100)==0
        display(Single);
    end;
    state=2;
    cycle=0;
    preState=state;
    flag1=1;
    while flag1==1
        dt=AA(StartPos+cycle+1,2)-AA(StartPos+cycle,2);
        switch state
            case 1
                Kt2=KCalculation(1,2,AA(StartPos+cycle,1));
                Kt4=KCalculation(1,4,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt2)
                    state=2;
                end;
                if rand(1)<(dt*Kt4)
                    state=4;
                end;
            case 2
                Kt3=KCalculation(2,3,AA(StartPos+cycle,1));
                Kt1=KCalculation(2,1,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt3)
                    state=3;
                end;
                if rand(1)<(dt*Kt1)
                    state=1;
                end;
            case 3
                Kt4=KCalculation(3,4,AA(StartPos+cycle,1));
                Kt2=KCalculation(3,2,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt4)
                    state=4;
                end;
                if rand(1)<(dt*Kt2)
                    state=2;
                end;
            case 4
                Kt1=KCalculation(4,1,AA(StartPos+cycle,1));
                Kt3=KCalculation(4,3,AA(StartPos+cycle,1));
                if rand(1)<(dt*Kt1)
                    state=1;
                end;
                if rand(1)<(dt*Kt3)
                    state=3;
                end;
        end;
%         if state~=preState
%             state
%         end;
        if state==4
            if preState==3
                Count(1)=Count(1)+1;
                TimeDisP(Single)=AA(StartPos+cycle,2);
            else
                TimeDisN(Single)=AA(StartPos+cycle,2);
                Count(2)=Count(2)+1;
            end;
            flag1=0;
        end;
        preState=state;
        cycle=cycle+1;
        if cycle>(length(AA)-StartPos-10)
            flag1=0;
        end;
    end;
end;

figure; bar(Count);
axis([0 3 0 ExpNum]);
xlabel('clockwise or anti-clockwise');
ylabel('Number of cases');

TimeDisP(find(TimeDisP==0))=[];
TimeDisP=TimeDisP-26;
figure; hist(TimeDisP,40);
xlabel('Dwell Time / h');
ylabel('Number of cases');

TimeDisN(find(TimeDisN==0))=[];
TimeDisN=TimeDisN-26;
figure; hist(TimeDisN,40);
xlabel('Dwell Time / h');
ylabel('Number of cases');
ylabel('Number of cases');
